package com.zyf.EasyNet.util.NeuralNetwork;

import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import java.text.DecimalFormat;
import java.util.Arrays;
import java.util.Random;

/**
 * @first_author zyflzz
 * @gmt_created 2022/5/31
 * @gmt_modified 2022/6/1
 */
public class BpNetUtils {

    /*
    约定：
        k 表示传入时样本序号     k = 0 1 2 ... K
        n 表示学习时迭代序号     n = 0 1 2 ... N
        i 表示输出层结点序号     i = 0 1 2 ... I
        h 表示隐含层结点序号     h = 0 1 2 ... H
        o 表示输出层结点序号     o = 0 1 2 ... O

        K 表示样本个数
        N 表示最大迭代数
        I 表示输入层结点个数
        H 表示隐含层节点个数
        O 表示输出层节点个数

        注：作为私有变量时有其他容易识别的名字，作为方法内部变量时按照这个规则来
     */

    /*----------------------------------------------------------------------------------------------------------*/

    /* 工具变量 */
    // 日志
    private static final Logger LOG = LoggerFactory.getLogger(BpNetUtils.class);
    // 设置保留小数点后6位
    DecimalFormat df = new DecimalFormat("0.000000");
    // 数组操作工具
    private static ArrayUtils arrayUtils = new ArrayUtils();
    // 激活函数库
    private static ActFuncUtils actFuncUtils = new ActFuncUtils();
    /*----------------------------------------------------------------------------------------------------------*/

    /* 基本输入变量 */
    // 样本矩阵 K行I列 下标ki，需要标准化后输入
    private double[][] sample;
    // 期望矩阵 K行O列 下标ko，需要标准化后输入
    private double[][] except;
    /*----------------------------------------------------------------------------------------------------------*/

    /* 学习调控变量 */
    // 期望误差
    private double epsilon;
    // 最大迭代次数
    private int MAX;
    // 学习因子 mu
    private double mu;
    // 学习因子 eta
    private double eta;
    /*----------------------------------------------------------------------------------------------------------*/

    /* 基本信息记录 */
    // 样本数
    private int sampleN;
    // 输入层结点数
    private int inputN;
    // 隐含层结点个数
    private int hideN;
    // 输出层结点数
    private int outputN;
    // 当前全局误差
    private double Error;
    // 当前迭代次数
    private int iterN;
    /*----------------------------------------------------------------------------------------------------------*/

    /* 权重矩阵 */
    // 隐含层到输入层之间的权重 p行n列 下标hi
    private double[][] thetaHI;
    // 输入层到隐含层之间的权重 n行p列 下标ih
    private double[][] thetaIH;
    // 输出层到隐含层之间的权重 m行p列 下标oh
    private double[][] thetaOH;
    // 隐含层到输出层之间的权重 p行m列 下标ho
    private double[][] thetaHO;
    /*----------------------------------------------------------------------------------------------------------*/

    /* 阈值矩阵 */
    // 隐含层阈值 p列 下标h
    private double[] thresholdH;
    // 输出层阈值 m列 下标o
    private double[] thresholdO;
    /*----------------------------------------------------------------------------------------------------------*/

    /* 中间迭代矩阵 */
    // 输入层节点的值 K行I列 下标i
    private double[][] inputLayer;
    // 隐含层节点的输入值 K行H列 下标kh
    private double[][] hideLayerI;
    // 隐含层节点的输出值 K行H列 下标kh
    private double[][] hideLayerO;
    // 输出层节点的输入值 K行O列 下标ko
    private double[][] outputLayerI;
    // 输出层节点的输出值 K行O列 下标ko
    private double[][] outputLayerO;
    /*----------------------------------------------------------------------------------------------------------*/

    /* 中间计算矩阵 */
    // 误差函数对输出层的偏导数 k行m列 下标o
    private double[][] deltaO;
    // 误差函数对隐含层的偏导数 k行p列 下标h
    private double[][] deltaH;
    /*----------------------------------------------------------------------------------------------------------*/

    /* 辅助变量 */
    // 用于反序列化
    private double[] sampleMax;
    private double[] sampleMin;
    private double[] sampleRange;
    private double[] exceptMax;
    private double[] exceptMin;
    private double[] exceptRange;
    /*----------------------------------------------------------------------------------------------------------*/

    /* 数据处理部分 */

    /**
     * 数据归一化
     * 传入一个需要归一化的矩阵，矩阵的每一行位一个样本，每一列位样本的变量
     * 归一化操作需要对同一类变量操作，也就是每一列变量操作
     * 不仅要传入一个归一化矩阵，还需要传入三个向量用来存储每一类变量的最值和极差，用于最后逆归一化
     * 归一化后会产生一个新的矩阵，以保证传入的数据不被破坏
     *
     * @param data       数据矩阵
     * @param maxArray   列最大值向量
     * @param minArray   列最小值向量
     * @param rangeArray 列极差向量
     */
    private double[][] norData(double[][] data, double[] maxArray, double[] minArray, double[] rangeArray) {
        // 首先要将矩阵转置，从而归类每一类变量
        double[][] dataT = arrayUtils.transpose(data);
        // 临时变量，用于存储最值
        double max, min;
        // 获取变量个数，即每一个样本有多少种变量，亦即数据矩阵的列数
        int c = dataT.length;
        // 遍历每一类变量，先存储好最值和极差向量
        for (int i = 0; i < c; i++) {
            // 首先获取最值
            max = arrayUtils.findMaxOrMin(dataT[i], true);
            min = arrayUtils.findMaxOrMin(dataT[i], false);
            // 设置最值向量和极差向量
            maxArray[i] = max;
            minArray[i] = min;
            rangeArray[i] = max - min;
        }
        // 调用归一化函数，进行归一化，其参数是数据矩阵的转置矩阵
        double[][] newData = arrayUtils.normalizeMatrix(dataT);
        // 归一化结束后需要再次转置过来
        return arrayUtils.transpose(newData);
    }

    /**
     * 数据逆归一化
     * 在学习完毕之后的模型运用时，传入的数据按照之前的规则进行归一化，输出的数据按照逆归一化解析成我们需要的数据
     * 在逆归一化时，需要根据输入的最小值向量和极差向量来逆归一化
     * 逆归一化后会产生一个新的矩阵，以保证传入的数据不被破坏
     *
     * @param data  输出数据矩阵
     * @param min   最小值向量
     * @param range 极差向量
     * @return 逆归一化后的数据
     */
    private double[][] invNorData(double[][] data, double[] min, double[] range) {
        // 获取行数，即样本数
        int k = data.length;
        // 获取列数，即变量数
        int N = data[0].length;
        // 新创建一个矩阵，用于存储逆归一化数据
        double[][] newData = new double[k][N];
        for (int i = 0; i < k; i++) {
            for (int j = 0; j < N; j++) {
                newData[i][j] = data[i][j] * range[j] + min[j];
            }
        }
        return newData;
    }
    /*----------------------------------------------------------------------------------------------------------*/

    /* 构造函数部分 */

    /**
     * 构造函数，需要传入样本矩阵、输入层的结点数和期望矩阵
     *
     * @param sample 样本矩阵
     * @param H      输入层结点数
     * @param except 期望矩阵
     */
    public BpNetUtils(double[][] sample, int H, double[][] except) {
        LOG.info("BP neural network simulation preparation");

        // 传入的参数不能位空
        if (sample == null || except == null) {
            LOG.error("sample or except matrix is null");
        } else {
            // 获取结点数和样本数
            int K = sample.length;
            int I = sample[0].length;
            int O = except[0].length;

            /* 输入和输出的初始化必须在这里先完成 */
            // 设置节点结点数和样本数
            this.sampleN = K;
            this.inputN = I;
            this.hideN = H;
            this.outputN = O;

            // 初始化逆归一化矩阵
            this.sampleMax = new double[I];
            this.sampleMin = new double[I];
            this.sampleRange = new double[I];
            this.exceptMax = new double[O];
            this.exceptMin = new double[O];
            this.exceptRange = new double[O];

            // 设置系统样本和期望矩阵
            this.sample = norData(sample, sampleMax, sampleMin, sampleRange);
            this.except = norData(except, exceptMax, exceptMin, exceptRange);

            // 设置系统输入矩阵并初始化
            this.inputLayer = new double[K][I];
            initInput(this.inputLayer, this.sample);

            // 根据样本数、输入层结点数、隐含层结点数、输出层结点数进行其他变量的基本初始化

            // 初始化隐含层输入矩阵和输出矩阵
            this.hideLayerI = new double[K][H];
            this.hideLayerO = new double[K][H];
            // 初始化输出层输入矩阵和输出矩阵
            this.outputLayerI = new double[K][O];
            this.outputLayerO = new double[K][O];
            // 初始化误差函数对隐含层和输出层的偏导数矩阵
            this.deltaO = new double[K][O];
            this.deltaH = new double[K][H];

            // 初始化权重矩阵
            this.thetaHI = new double[H][I];
            initTheta(thetaHI);
            this.thetaIH = new double[I][H];
            thetaIH = arrayUtils.transpose(thetaHI);
            this.thetaHO = new double[H][O];
            initTheta(thetaHO);
            this.thetaOH = new double[O][H];
            thetaHO = arrayUtils.transpose(thetaOH);
            // 初始化阈值矩阵
            this.thresholdH = new double[H];
            Arrays.fill(thresholdH, 0.0);
            this.thresholdO = new double[O];
            Arrays.fill(thresholdO, 0.0);

            // 初始化精度、迭代数、学习效率
            this.epsilon = 0.03;
            this.MAX = 1000;
            this.eta = 0.05;
            this.mu = 0.05;

            LOG.info("BP neural network simulation preparation is OK");
        }
    }
    /*----------------------------------------------------------------------------------------------------------*/

    /* 主要功能模块 */

    /**
     * 学习迭代
     * 返回值为 0~9 数字组成的代码串
     * 具体返回值      含义
     * 0       执行成功，以误差精度满足条件迭代结束
     * 1       执行成功，以迭代次数满足条件迭代结束
     * 2       执行失败，样本数组为空
     * 3       执行失败，输入层-隐含层权重为空
     * 4       执行失败，隐含层-输出层权重为空
     * 5       执行失败，隐含层阈值为空
     * 6       执行失败，输出层阈值为空
     *
     * @return 执行情况
     */
    public String bpIteration() {
        LOG.info("BP neural network simulation start iteration");
        // 检查数据是否不完整
        String checkResult = check();
        if (!checkResult.equals("")) {
            return checkResult;
        }
        // 记录输出代码
        String code = "";
        // 获取样本个数、输入层结点数、隐含层结点数、输出层结点数、最大迭代数
        int I = inputN;
        int H = hideN;
        int O = outputN;
        int K = sampleN;
        int N = MAX;
        // 用于计算全局误差
        double E_all;
        // 用于计算误差
        double E;
        // 用于记录迭代次数
        int n;
        // 第一层迭代 学习轮数
        for (n = 0; n < N; n++) {
            // 每一轮学习都要更新全局误差
            E_all = 0.0;
            // 第二层迭代 样本选择
            for (int k = 0; k < K; k++) {
                /* 对于选定的样本，第一步需要计算各层的输入输出 */
                // 首先是隐含层的输入和输出
                for (int h =0; h < H; h++) {
                    // 根据该隐含层结点对于输入层各结点的权重、输入层各节点的值、该隐含层结点的阈值来计算该隐含层结点的输入值
                    hideLayerI[k][h] = calcInput(thetaHI[h], inputLayer[k], thresholdH[h]);
                    // 根据该隐含层结点的输入值计算该隐含层结点的输出值
                    hideLayerO[k][h] = calcOutput(hideLayerI[k][h]);
                }
                // 其次是输出层的输入和输出
                for (int o = 0; o < O; o++) {
                    // 根据该输出层结点对于隐含层各结点的权重、隐含层各结点的值、该输出层结点的阈值来计算该输出层结点的输入值
                    outputLayerI[k][o] = calcInput(thetaOH[o], hideLayerO[k], thresholdO[o]);
                    // 根据该输出层结点的输入值计算该输出层结点的输出值
                    outputLayerO[k][o] = calcOutput(outputLayerI[k][o]);
                }

                /* 对于已经计算好隐含层和输出层的输入和输出，第二步需要计算误差函数对于隐含层和输出层的偏导数，即梯度 */
                // 首先计算对输出层的偏导数，即输出层的梯度
                for (int o = 0; o < O; o++) {
                    // 根据该输出层的输出和期望值来激素那输出层梯度
                    deltaO[k][o] = deltaO(outputLayerO[k][o], except[k][o]);
                }
                // 其次计算对隐含层的偏导数，即隐含层的梯度
                for (int h = 0; h < H; h++) {
                    // 根据该隐含层的输入、该隐含层到各输出层的权重和各输出层的梯度来计算隐含层的梯度
                    deltaH[k][h] = deltaH(hideLayerI[k][h], thetaHO[h], deltaO[k]);
                }

                /* 对于已经计算好的梯度，就可以利用梯度下降法修正权值 */
                // 首先新创建两个矩阵，用来存储每次修正的增量
                double[][] DeltaHO = new double[H][O];
                double[][] DeltaIH = new double[I][H];
                // 先是隐含层和输出层之间的权重
                for (int o = 0; o < O; o++) {
                    for (int h = 0; h < H; h++) {
                        DeltaHO[h][o] = mu * deltaO[k][o] * hideLayerO[k][h];
                        thetaHO[h][o] += DeltaHO[h][o];
                        thetaOH[o][h] += DeltaHO[h][o];
                    }
                }
                // 然后是输入层和隐含层之间的权重
                for (int h = 0; h < H; h++) {
                    for (int i = 0; i < I; i++) {
                        DeltaIH[i][h] = eta * deltaH[k][h] * inputLayer[k][i];
                        thetaIH[i][h] += DeltaIH[i][h];
                        thetaHI[h][i] += DeltaIH[i][h];
                    }
                }
                /* 最后计算这一个样本的误差 */
                E = E(except[k], outputLayerO[k]);
                E_all += E;
            }
            // 每一轮结束后，所有样本都已经遍历，需要计算全局误差
            E_all /= K;
            this.Error = E_all;
            if (E_all < epsilon) {
                code = "0";
                n++;
                break;
            }
        }
        this.iterN = n;
        if (n >= MAX) {
            code = "1";
        }
        LOG.info("End simulation of BP neural network");
        return code;
    }

    /**
     * 对输入数据归一化，注意这个和训练数据不一样
     *
     * @param input 输入数据
     * @return 归一化后的数据
     */
    private double[][] norInput(double[][] input) {
        // 输入的 k 组数据
        int k = input.length;
        // 每组数据的变量个数
        int n = except[0].length;
        double[][] newInput = new double[k][n];
        for (int i = 0; i < k; i++) {
            for (int j = 0; j < n; j++) {
                newInput[i][j] = (input[i][j] - sampleMin[j]) / sampleRange[j];
            }
        }
        return newInput;
    }

    /**
     * 应用模型
     *
     * @param in 输入数据
     * @return 输出结果
     */
    public double[][] applyModel(double[][] in) {
        // 输入的 k 组数据
        int k = in.length;
        // 输出向量个数
        int m = except[0].length;
        // 将输入归一化
        double[][] input = norInput(in);
        double hideI;
        double outI;
        double[][] output = new double[k][m];
        double[][] hideLayer = new double[k][hideN];
        for (int i = 0; i < k; i++) {
            for (int p = 0; p < hideN; p++) {
                hideI = 0.0;
                for (int j = 0; j < input[0].length; j++) {
                    hideI += input[i][j] * thetaHI[p][j];
                }
                hideLayer[i][p] = calcOutput(hideI);
            }
            for (int o = 0; o < m; o++) {
                outI = 0.0;
                for (int p = 0; p < hideN; p++) {
                    outI += hideLayer[i][p] * thetaOH[o][p];
                }
                output[i][o] = calcOutput(outI);
            }
        }
        return invNorData(output, exceptMin, exceptRange);
    }
    /*----------------------------------------------------------------------------------------------------------*/

    /* 核心计算函数模块 */

    /**
     * 计算隐含层/输出层单个节点的输入值
     *
     * @param theta     输入层到隐含层该节点的权值/隐含层到输出层该节点的权值
     * @param layer     输入层的各个节点/隐含层的各个节点
     * @param threshold 隐含层的阈值/输出层的阈值
     * @return 输入值
     */
    private double calcInput(double[] theta, double[] layer, double threshold) {
        int len = layer.length;
        double result = 0.0;
        for (int l = 0; l < len; l++) {
            result += theta[l] * layer[l];
        }
        result -= threshold;
        return result;
    }

    /**
     * 计算隐含层/输出层单个节点的输出值
     *
     * @param param 隐含层/输出层的输入值
     * @return 输出值
     * <p>
     * {@link BpNetUtils#f(double)}
     */
    private double calcOutput(double param) {
        return f(param);
    }

    /**
     * 计算误差函数对输出端输入的偏导数
     *
     * @param inputO 输出端输入
     * @param except 输出端阈值
     * @return 误差函数对输出端输入的偏导数
     * <p>
     * {@link BpNetUtils#f(double)}
     * {@link BpNetUtils#df(double)}
     */
    private double deltaO(double inputO, double except) {
        return (except - f(inputO)) * df(inputO);
    }

    /**
     * 计算误差函数对隐含层输入的偏导数
     *
     * @param inputH 隐含层输入
     * @param thetaO 隐含层到输出层该节点的权重
     * @param deltaO 误差函数对输出端输入的偏导数
     * @return 误差函数对隐含层输入的偏导数
     * <p>
     * {@link BpNetUtils#deltaO(double, double)}
     * {@link BpNetUtils#df(double)}
     */
    private double deltaH(double inputH, double[] thetaO, double[] deltaO) {
        int m = thetaO.length;
        double result = 0.0;
        for (int o = 0; o < thetaO.length; o++) {
            result += thetaO[o] * deltaO[o];
        }
        result *= df(inputH);
        return result;
    }

    /**
     * 计算输出误差
     *
     * @param except      期望值
     * @param outputLayer 实际输出
     * @return 平方误差
     */
    private double E(double[] except, double[] outputLayer) {
        int m = except.length;
        double err;
        double sum = 0.0;
        for (int o = 0; o < m; o++) {
            err = Math.abs(except[o] - outputLayer[o]);
            sum += err * err;
        }
        return 0.5 * sum;
    }

    /**
     * 激活函数： tanh 函数
     *
     * @param param 函数参数
     * @return 函数值
     */
    private double f(double param) {
        return actFuncUtils.tanh(param);
    }

    /**
     * 激活函数导数
     *
     * @param param 函数参数
     * @return 函数值
     */
    private double df(double param) {
        return actFuncUtils.dTanh(param);
    }
    /*----------------------------------------------------------------------------------------------------------*/

    /* 辅助函数模块 */

    /**
     * 初始化权重矩阵
     *
     * @param theta 权重矩阵
     */
    private void initTheta(double[][] theta) {
        Random random = new Random();
        for (int i = 0; i < theta.length; i++){
            for (int j = 0; j < theta[0].length; j++){
                theta[i][j] = random.nextDouble();
            }
        }
    }

    /**
     * 初始化输入矩阵
     *
     * @param input  输入矩阵
     * @param sample 样本矩阵
     */
    private void initInput(double[][] input, double[][] sample) {
        for (int i = 0; i < sample.length; i++) {
            System.arraycopy(sample[i], 0, input[i], 0, sample[0].length);
        }
    }

    /**
     * 数据检查
     *
     * @return 检查结果代码
     */
    private String check() {
        // 返回代码
        StringBuilder code = new StringBuilder("");
        if (sample == null) {
            code.append("2");
        }
        if (thetaHI == null || thetaIH == null) {
            code.append("3");
        }
        if (thetaHO == null || thetaOH == null) {
            code.append("4");
        }
        if (thresholdH == null) {
            code.append("5");
        }
        if (thresholdO == null) {
            code.append("6");
        }
        return String.valueOf(code);
    }
    /*----------------------------------------------------------------------------------------------------------*/

    /* getter/setter 类 */

    public double[][] getSample() {
        return sample;
    }

    public void setSample(double[][] sample) {
        arrayUtils.normalizeMatrix(sample);
        this.sample = sample;
    }

    public double[][] getExcept() {
        return except;
    }

    public void setExcept(double[][] except) {
        this.except = except;
    }

    public double getEpsilon() {
        return epsilon;
    }

    public void setEpsilon(double epsilon) {
        this.epsilon = epsilon;
    }

    public int getMAX() {
        return MAX;
    }

    public void setMAX(int MAX) {
        this.MAX = MAX;
    }

    public int getHideN() {
        return hideN;
    }

    public void setHideN(int hideN) {
        this.hideN = hideN;
    }

    public double[][] getThetaHI() {
        return thetaHI;
    }

    public void setThetaHI(double[][] thetaHI) {
        this.thetaHI = thetaHI;
        this.thetaIH = arrayUtils.transpose(thetaHI);
    }

    public double[][] getThetaIH() {
        return thetaIH;
    }

    public void setThetaIH(double[][] thetaIH) {
        this.thetaIH = thetaIH;
        this.thetaHI = arrayUtils.transpose(thetaIH);
    }

    public double[][] getThetaOH() {
        return thetaOH;
    }

    public void setThetaOH(double[][] thetaOH) {
        this.thetaOH = thetaOH;
        this.thetaHO = arrayUtils.transpose(thetaOH);
    }

    public double[][] getThetaHO() {
        return thetaHO;
    }

    public void setThetaHO(double[][] thetaHO) {
        this.thetaHO = thetaHO;
        this.thetaOH = arrayUtils.transpose(thetaHO);
    }

    public double[] getThresholdH() {
        return thresholdH;
    }

    public void setThresholdH(double[] thresholdH) {
        this.thresholdH = thresholdH;
    }

    public double[] getThresholdO() {
        return thresholdO;
    }

    public void setThresholdO(double[] thresholdO) {
        this.thresholdO = thresholdO;
    }

    public double getMu() {
        return mu;
    }

    public void setMu(double mu) {
        this.mu = mu;
    }

    public double getEta() {
        return eta;
    }

    public void setEta(double eta) {
        this.eta = eta;
    }
    /*----------------------------------------------------------------------------------------------------------*/

    /* getter 类 */
    public double getError() {
        return Error;
    }

    public int getIterN() {
        return iterN;
    }

    public double[][] getOutputLayerO() {
        return outputLayerO;
    }

    public double[] getSampleMax() {
        return sampleMax;
    }

    public double[] getSampleMin() {
        return sampleMin;
    }

    public double[] getSampleRange() {
        return sampleRange;
    }

    public double[] getExceptMax() {
        return exceptMax;
    }

    public double[] getExceptMin() {
        return exceptMin;
    }

    public double[] getExceptRange() {
        return exceptRange;
    }
}